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Introduction 


This  document  describes  a  set  of  signal-processing  algorithms,  as  implemented  by  the  Computer 
Engineering  Research  Laboratory  at  Georgia  Tech.  The  routines  are  presented  as  a  representative 
collection  of  operations  for  processing  Infrared  Focal-Plane  Array  signals. 

For  the  purposes  of  testing  and  dissemination,  each  algorithm  is  presented  as  a  stand-alone 
FORTRAN  program.  These  programs  are  based  upon  a  core  harness  routine  which  supports  the 
input/output  of  a  common  data  format  (Georgia  Tech  Algorithm  Evaluation  Data  Format  - 
described  in  the  Harness  section).  The  modular  implementations  offer  several  benefits: 

0  simplification  of  the  generation  of  test  vectors  for  the  verification  of  alternate 
implementations 

°  capability  for  testing  various  algorithm  combinations,  without  re-compilation 

0  support  for  multiple  language  and/or  processor-platform  implementations 

II.  Harness 

A.  Description 

The  Harness  program  shown  below  is  the  basis  of  the  input/output  methodology  used  by  all  of 
the  routines  in  this  document.  The  code  implements  a  simple  Pass-Through  module  which  reads 
a  data  stream,  picking  off  the  FPA  pixel  data,  and  writing  the  data  onto  an  output  data  stream. 

The  Georgia  Tech  Algorithm  Evaluation  Data  format  is  a  simple  ASCII  text  representation  of  a 
data  stream.  The  data  stream  has  two  major  components  -  the  Field  Header  and  the  Field  Data. 
The  harness  of  each  module  processes  the  data  stream  by  reading  each  line  and  checking  for 
Field  Headers  which  are  relevant' to  that  module.  Any  lines  which  are  not  relevant,  or 
unrecognized,  are  immediately  placed  upon  the  output  data  stream.  As  soon  as  a  relevant  Field 
Header  is  recognized,  the  Field  Data  which  follows  is  processed  in  a  manner  which  is  appropriate 
to  that  module  and  Field  Header.  This  scheme  provides  for  the  chaining  of  modules  output-to- 
input,  without  either  module  requiring  knowledge  of  all,  or  any,  of  the  other  module's  data 
formats.  In  typical  use,  controls  for  many  modules  could  be  included  in  a  single  data  stream; 
each  module  would  only  process  data  intended  for  it.  For  example,  suppose  a  test  setup  was 
composed  of  the  following  pipeline: 

Input  data  stream  — >  Spatial  Filter  - — >  Simple  Threshold  - — >  Output  Data  Stream 


Harness 


The  data  stream  might  appear  as  follows: 


Input  Data  Stream 

Description 

Used  by 

Action 

Dimensions 

128 

Reid  Header 

Reid  Data 

Spatial  Filter  and 
Simple  Threshold 

input 

Simple  Thresholding  Limits 

0  256 

Field  Header 

Reid  Data 

Simple  Threshold 

input 

Pixel  Data 

99  93  ... 

81  76  ... 

Field  Header 

Field  Data 

Field  Data 

Spatial  Filter  and 
Simple  Threshold 

input,  modified 

End 

Reid  Header 

Spatial  Filter  and 
Simple  Threshold 

input 

Output  Data  Stream 

Description 

Generated  by 

Action 

Dimensions 

128 

Field  Header 

Field  Data 

Spatial  Filter  and 
Simple  Threshold 

copied  to  output 
data  stream 

Simple  Thresholding  Limits 

0  256 

Field  Header 

Field  Data 

Simple  Threshold 

copied  to  output 
data  stream 

Simple  Thresholding 

Statistics 

0  256  1024 

Reid  Header 

Reid  Data 

Simple  Threshold 

generated,  then 
placed  on  output 
data  stream 

Pixel  Data 

99  93  ... 

81  76  ... 

Field  Header 

Reid  Data 

Reid  Data 

Spatial  Filter  and 
Simple  Threshold 

modified,  then 
placed  on  output 
data  stream 

End 

Field  Header 

Spatial  Filter  and 
Simple  Threshold 

copied  to  output 
data  stream 

B.  Module  Listing 

PROGRAM  HARNESS 
C 

C  This  program  acts  as  a  harness  for  testing  various 

C  Signal  Processing  Routines 

C 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 
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Harness 


400  Tenth  St.  CRB  390 
Atlanta,  GA  30332-0540 
Contact:  Andrew  Henshaw  (404)894-2521 

Written  by  A.  M.  Henshaw  Jan  23,  1990 
Using  Microsoft  Fortran 

CHARACTER* (*)  Dim,  Pixels  . 

PARAMETER  (Dim  “'Dimensions') 

PARAMETER  (Pixels* 'Pixel  Data') 

PARAMETER  (maxSize-64) 

INTEGER  n 

INTEGER  in (maxSize,maxSize)  ,  out(maxSize,  maxSize) 

CHARACTER  header *7 2 
LOGICAL  runFlag 

WRITE  (6,*)  '%  Processed  by  Pass  Thru  module.' 
runFlag  =  . TRUE . 

DO  WHILE  (runFlag) 

READ  (5,1000)  header 
1000  FORMAT  (A72) 

IF  (header. EQ. Dim)  THEN 
READ  (5,*)  n 
WRITE  (6,*)  Dim 
WRITE  (6,*)  n 

ELSE  IF  (header. EQ. Pixels)  THEN 

READ  (5,*)  ( (in (row, col) , col=l, n) , row=l, n) 

CALL  PassThru  (n,  in,  out) 

WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (out (row, col) , col=l, n) , row=l, n) 

ELSE  IF  ( header. EQ. 'End' )  THEN 
WRITE  (6,*)  'End' 
runFlag  =  .FALSE. 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 

END 

C* ******************************************************************* 

SUBROUTINE  PassThru  (n,  in,  out) 

PARAMETER  (maxSize=64) 

INTEGER  n,  row,  col 
INTEGER  in (maxSize,  maxSize) 

INTEGER  out (maxSize,  maxSize) 


DO  30  row  =  1,  n 
DO  30  col  =  1,  n 

out (row, col)  =  in (row, col) 
30  CONTINUE 

RETURN 

END 
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Non-Uniformity  Compensation 


m.  Non-Uniformitv  Compensation 

A.  Description 

The  non-uniformity  compensation  algorithm  provides  a  pixel-by-pixel  correction  of  the  actual 
pixel  response  to  the  desired  response.  The  current  algorithm  uses  up  to  a  five-point,  piecewise- 
linear  correction  to  the  pixel  intensity.  The  correction  is  determined  by  sending  a  specified 
number  of  calibration  frames  through  the  process.  Each  of  the  calibration  frames  will  have  been 
generated  by  exposing  the  focal-plane  array  to  a  known  intensity  so  that  a  desired  pixel  intensity 
is  expected  at  each  pixel. 

Dead,  or  inadequately  responsive,  pixels  are  assumed  to  have  been  marked  by  another  module 
and,  during  processing,  they  are  replaced  by  the  intensity  of  the  previous  pixel. 

After  calibration  is  performed,  the  algorithm  enters  the  processing  phase.  For  each  pixel  which  is 
to  be  processed,  it  is  first  determined  if  it  is  a  dead  pixel.  For  normal  pixels,  the  calibration 
intensities  are  searched  to  determine  which  section  should  be  used  for  correction.  After  the 
section  is  determined,  a  linear  interpolation  is  performed  using  the  input  pixel  intensity  to 
interpolate  between  the  desired  responses. 

B.  Data  fields 


Action 

Field  Header 

Field  Data 

Data  Type 

input 

Dimensions 

FPA  dimension 

Integer 

input 

Calibration  Frames 

Count  of  calibration 
frames 

Integer 

input 

Calibration  Input 

Integer 
fl..  Count] 

input 

Calibration  Pixel  Data 

Array  of  pixel 
response  data  for  one 
input  reference 

Integer 

[1..  Dimension] 

[1.. Dimension] 

modify 

Pixel  Data 

Pixel  data  array 

Integer 

[1.. Dimension] 
fl.. Dimension! 

C.  Module  Listing 

PROGRAM  NUNICOMP 
C 

C  Non-Uniform  Compensation  Test  Module 

C 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 

C  400  Tenth  St.  CRB  390 

C  Atlanta,  GA  30332-0540 

C  Contact:  Andrew  Henshaw  (404)894-2521 

C 

C  conforms  to  the  Ga .  Tech  Algorithm  Evaluation  Data  Format 
C 

C  Fortran  translation  of  Occam  code 

C  Steve  Gieseking 
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Non-Uniformity  Compensation 


Roy  W.  Melton  Feb  1,  1990 

Harness  written  by  Andrew  Henshaw  Jan  23,  1990 
Using  Microsoft  Fortran 

CHARACTER* (* )  Callnp,  CalOutp,  Dim,  Pixels,  Sect 
Valid  Section  Headers 


C 


PARAMETER  (Callnp 
PARAMETER  (CalOutp 
PARAMETER  (Dim 
PARAMETER  (Pixels 
PARAMETER  (Sect 


'Calibration  Input') 
'Calibration  Pixel  Data') 
'Dimensions ' ) 

' Pixel  Data  1 ) 

'Calibration  Frames') 


PARAMETER  (maxSize  =  128) 
PARAMETER  (maxCalFrames  =  5) 


!  maximum  FPA  size 
!  default  value 


INTEGER  N,  Count,  Sections 
INTEGER  Ic  (maxCalFrames) 

INTEGER  In  (maxSize,  maxSize),  Out (maxSize,  maxSize) 
INTEGER  Oc  (maxCalFrames,  maxSize,  maxSize) 
CHARACTER  Header *72 
LOGICAL  runFlag 


Count  *  1 

Sections  -  maxCalFrames 

WRITE  (6,*)  '%  Processed  by  Non-Uniformity  Compensation  Module.' 

runFlag  =  . TRUE . 

DO  WHILE  (runFlag) 

READ  (5,1000)  Header 
1000  FORMAT  (A72) 

IF  (Header. EQ. Callnp)  THEN 
IF  (Count. LE. Sections)  THEN 
READ  (5,*)  Ic  (Count) 

WRITE  (6,*)  Callnp 
WRITE  (6,*)  Ic  (Count) 

ELSE 

WRITE  (6,*)  Callnp 
END  IF 


ELSEIF  (Header. EQ. CalOutp)  THEN 
IF  (Count. LE. Sections)  THEN 

READ  (5,  *)  ((Oc  (Count,  Row,  Col),  Col=l,N) ,  Row-1, N) 
WRITE  (6,*)  CalOutp 

WRITE  (6,*)  ( (Oc  (Count,  Row,  Col),  Col=l,N),  Row-1, N) 

Count  =  Count  +  1 
ELSE 

WRITE  (6,*)  CalOutp 
ENDIF 

ELSEIF  (Header. EQ. Dim)  THEN 
READ  (5,*)  N 
WRITE  (6,*)  Dim 
WRITE  (6,*)  N 

ELSE  IF  (Header. EQ. Pixels)  THEN 

IF  ( (Count. GT.l) .AND. (Sections. GT.l) )  THEN 
READ  (5,*)  ( (In (row, col) , col-1, n) , row=l, n) 

IF  (Count. LE. Sections)  THEN 
Sections  =  Count  -  1 
ENDIF 

CALL  NonUniformityCompensation 
+  (In,  Out,  Oc,  Ic,  N,  Sections) 
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WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (Out (row, col) ,  col=l,n),  row=l,n) 

ELSE 

WRITE  (6,*)  Pixels 
ENDIF 

ELSEIF  (Header. EQ. Sect)  THEN 
READ  (5,*)  Sections 
WRITE  (6,*)  Sect 
WRITE  (6,*)  Sections 

ELSE  IF  (Header. EQ. 'End' )  THEN 
WRITE  (6,*)  'End' 
runFlag  =  .FALSE. 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 


END 


SUBROUT INE  NonUni f o rmit yCompe  n s  a t i on 
+  (In,  Out,  Oc,  Ic,  N,  Sections) 

PARAMETER  (MAXCALFRAMES=5 ) 

PARAMETER  (MAXSIZE=64) 

INTEGER  In  (MAXSIZE,  MAXSIZE) ,  Out  (MAXSIZE,  MAXSIZE) 
INTEGER  Oc  (MAXCALFRAMES,  MAXSIZE,  MAXSIZE) 

INTEGER  IC  (MAXCALFRAMES) 

INTEGER  N,  Sections 

INTEGER  I,  J,  LastPixel,  Section 


LastPixel  =  0 
DO  20  I  =  1,  N 
DO  20  J  =  1,  N 

IF  (Oc  (1,  I,  J)  .EQ.  65535)  THEN 
Out  (I,  J)  =  LastPixel 
ELSE 

Section  =  1 

10  IF  ( (Section. LT. (Sections  -  1)).AND. 

+  (In  (I,  J) .GE .Oc  ((Section  +1),  I,  J) ) )  THEN 

Section  =  Section  +  1 
GOTO  10 
ENDIF 


IF  (In  (I,  J)  . LT . Oc  (Section,  I,  J)  )  THEN 


Out 

(I,  J)  =  Oc 

(Section,  I,  J) 

ELSE 

Out 

(I,  J)  =  (In 

(I,  J)  -  OC  (Section, 

I,  J))  * 

(Ic 

(Section  +  1)  - 

Ic 

(Section) )  / 

(Oc 

( (Section  +  1 ) , 

I, 

J) 

- 

Oc 

(Section,  I,  J) 

)  + 

Ic 

(Section) 

ENDIF 

IF  (Out  (I,  J).GT. 

65535)  THEN 

Out  (I,  J)  =  65535 
ENDIF 

LastPixel  =  Out  (I,  J) 
ENDIF 

20  CONTINUE 
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RETURN 

END 


I. 

! 
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Spatial  Filtering 


EL _ Spatial  Filtering 

A.  Description 

The  spatial  filtering  algorithm  performs  a  convolution  of  the  image  with  a  3x3  coefficient  mask. 
This  implementation  supports  a  four  point  symmetric  mask.  Separate  masks  are  used  for  the 
edge  pixels  since  not  all  of  the  pixels  which  are  needed  are  defined.  This  allows  more  general 
application  of  boundary  conditions  than  would  be  available  if  the  undefined  pixels  were  treated 
as  zeros  and  the  same  mask  was  used. 

Since  the  filter  allows  negative  coefficients  in  the  mask,  it  is  possible  to  generate  negative  output 
intensities.  The  coding  allows  the  intensity  to  be  output  limited  to  a  positive  range. 

B.  Data  Fields 


Action 

Field  Header 

Field  Data 

Data  Type 

Dimensions 

FPA  dimension 

Integer 

input 

Spatial  Filter  Controls 

Filter  coefficients 
(Comer  coefficients) 

Integer 

n..4i 

Filter  coefficients 
(Top  coefficients) 

Integer 

[1-4] 

Integer 

[1-4] 

Filter  coefficients 
(Center  coefficients) 

Integer 

[1-4] 

modify 

Pixel  Data 

Pixel  data  array 

Integer 

[1..  Dimension] 

[1.. Dimension] 

C.  Module  Listing 

PROGRAM  SPFILT 
C 

C  Spatial  Filtering  Test  Module 

C 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 

C  400  Tenth  St.  CRB  390 

C  Atlanta,  GA  30332-0540 

C  Contact:  Andrew  Henshaw  (404)894-2521 

C 

C  conforms  to  the  Ga .  Tech  Algorithm  Evaluation  Data  Format 
C 

C  Fortran  translation  of  Occam  code 

C  Steve  Gieseking 

C  Roy  Melton 

C 

C  Harness  written  by  Andrew  Henshaw  Jan  23,  1990 

C  Using  Microsoft  Fortran 

C 

CHARACTER* (*)  Controls,  Dim,  Pixels 

PARAMETER  (Controls  =  'Spatial  Filter  Controls') 
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PARAMETER  (Dim  =  'Dimensions') 

PARAMETER  (Pixels  =  'Pixel  Data') 

PARAMETER  (maxSize  -  64) 

PARAMETER  (SF_CONTROL_SIZE  =  4) 

INTEGER  N 

INTEGER  In  (maxSize,  maxSize) ,  Out  (maxSize,  maxSize) 

INTEGER  C  ( SF_CONTROL_S I ZE ,  SF_CONTROL_SIZE) 

CHARACTER  header*72 
LOGICAL  runFlag 

WRITE  (6,*)  '%  Processed  by  Spatial  Filtering  Module.' 

CALL  DefaultFilterControls  (C) 
runFlag  =  . TRUE . 

DO  WHILE  (runFlag) 

READ  (5,1000)  header 
1000  FORMAT  (A72) 

IF  (header. EQ. Controls)  THEN 

READ  (5,  *)  ( (C  (row,  col),  col-1,  SF_CONTROL_SIZE) , 

+  row-1,  SF_CONTROL_SIZE  ) 

WRITE  (6,*)  Controls 

WRITE  (6,*)  ( (C  (row,  col),  col-1,  SF_CONTROL_SIZE) , 

+  row=l,  SF_CONTROL_S I ZE  ) 

ELSE  IF  (header. EQ. Dim)  THEN 
READ  (5,*)  N 
WRITE  (6,*)  Dim 
WRITE  (6,*)  N 

ELSE  IF  (header. EQ. Pixels)  THEN 

READ  (5,*)  ( (In (row, col) , col-1, n) , row=l, n) 

CALL  SpatialFilter  (In,  Out,  C,  N) 

WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (Out (row, col) , col-1, n) , row-1, n) 

ELSE  IF  ( header. EQ. 'End')  THEN 
WRITE  (6,*)  'End' 
runFlag  -  .FALSE. 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 

END 

C***Filter  Control  Defaults****************************************** 
SUBROUTINE  DefaultFilterControls  (Control) 

PARAMETER  (SF_CONTROL_SIZE  =  4) 

INTEGER  Control  (SF_CONTROL_SIZE,  SF_CONTROL_SIZE) 

INTEGER  I,  J 

DO  210  I  =  1,  4 
DO  200  J  =  1,  3 

Control  (I,  J)  =  0 
200  CONTINUE 

Control  (I,  4)  =  16384 
210  CONTINUE 

RETURN 

END 
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C***Spatial  Filter**************************** 
SUBROUTINE  SpatialFilter  (In,  Out,  C,  N) 

PARAMETER  (MAXSIZE  »  64) 

PARAMETER  (SF  CONTROL  SIZE  =  4) 


INTEGER  In  (MAXSIZE,  MAXSIZE),  Out  (MAXSIZE,  MAXSIZE) 
INTEGER  C  ( SF_CONTROL_S IZE ,  SF_CONTROL_SIZE) 

INTEGER  N 
INTEGER  I,  J 


DO  100  I  =  1,  N 
DO  100  J  =  1,  N 
IF  (I.EQ.l)  THEN 


IF  (J.EQ.l)  THEN 


Out 

(I,  J)  =  (In  (1+1,  J+l) 

*  C 

(lr 

1)  )  + 

+ 

(In  (I,  J+l) 

*  C 

(lr 

2))  + 

+ 

(In  <1+1,  J) 

*  c 

dr 

3))  + 

+ 

(In  (I,  J) 

*  c 

dr 

4)  ) 

ELSEIF 

(J.EQ.N)  THEN 

Out 

(I,  J)  -  (In  <1+1,  J-l) 

*  c 

dr 

D)  + 

+ 

(In  (I,  J-l) 

*  c 

dr 

2))  + 

+ 

(In  <1+1,  J) 

*  c 

dr 

3))  + 

+ 

(In  (I,  J) 

*  c 

dr 

4)  ) 

ELSE 

Out 

(I,  J)  =  ((In  <1+1,  J-l)  +  : 

In  (1+1,  J+l))  * 

C  (2 

,  1)  )  + 

+ 

((In  (I,  J-l) 

+  : 

In  (I 

,  J+l)) 

* 

C  (2 

,  2)  )  + 

+ 

(In  <1+1,  J) 

* 

C  (2 

,  3)  )  + 

+ 

(In  (I,  J) 

* 

C  (2 

,  4)) 

ENDIF 

ELSEIF  (I.EQ.N)  THEN 

IF  (J. 

EQ.l)  THEN 

Out 

(I,  J)  -  (In  <1-1,  J+l) 

*  c 

(lr 

D)  + 

+ 

(In  (I,  J+l) 

*  c 

(lr 

2))  + 

+ 

(In  <1-1,  J) 

*  c 

(If 

3))  + 

+ 

(In  (I,  J) 

*  c 

dr 

4)  ) 

ELSEIF 

(J.EQ.N)  THEN 

Out 

(I,  J)  -  (In  <1-1,  J-l) 

*  c 

(If 

D)  + 

+ 

(In  (I,  J-l) 

*  c 

(1, 

2))  + 

+ 

(In  (1-1,  J) 

*  c 

d. 

3))  + 

+ 

(In  (I,  J) 

*  c 

(1, 

4)  ) 

ELSE 

Out 

(I,  J)  =  ((In  (1-1,  J-l)  + 

In  (I 

-1,  J+l))  * 

C  (2 

,  1)  )  + 

+ 

((In  (I,  J-l) 

+ 

In  (I 

,  J+l)) 

★ 

C  (2 

,  2)  )  + 

+ 

(In  <1-1,  J) 

★ 

C  (2 

,  3)  )  + 

+ 

(In  (I,  J) 

★ 

C  (2 

,  4)) 

ENDIF 

ELSEIF  (J.EQ.l)  THEN 

Out  (I 

,  J)  =  ((In  (1-1,  J+l) 

+  In 

(1+1 

,  J+l)) 

*  C 

(3, 

D)  + 

+ 

(In  (I,  J+l) 

*  c 

(3, 

2)  )  + 

+ 

((In  <1-1,  J) 

+  In 

(I  +1,  J)) 

*  c 

(3, 

3)  )  + 

+ 

(In  (I,  J) 

*  c 

(3, 

4)  ) 

ELSEIF  (J.EQ.N)  THEN 

Out  (I 

,  J)  =  ((In  <1-1,  J-l) 

+  In 

(1+1 

,  J-l)) 

*  c 

(3, 

1)  )  + 

+ 

(In  (I,  J-l) 

*  c 

(3, 

2)  )  + 

+ 

((In  <1-1,  J)  + 

In  (1+1, 

J)  ) 

*  c 

(3, 

3)  )  + 

+ 

(In  (I,  J) 

*  c 

(3, 

4)  ) 

ELSE 

Out  (I 

,  J)  =  ((In  (1-1,  J-l) 

+  In 

(1+1 

J-l) 

+ 

+ 

In  (1-1,  J+l) 

+  In 

(1+1 

J+l) 

)  * 

C  (4, 

1)  )  + 

+ 

((In  (I,  J-l) 

+  In 

(1, 

J+l)  ) 

*  1 

C  (4, 

2)  )  + 

+ 

((In  (1-1,  J) 

+  In 

(1  +  1 

.,  J)  ) 

*  1 

C  (4, 

3)  )  + 

+ 

(In  (I,  J) 

*  1 

C  (4, 

4)  ) 

ENDIF 
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Out  (I,  J)  -  Out  (I,  J)  /  16384 
IF  (Out  (I,  J)  .  LT .  0 )  THEN 
Out  (I,  J)  -  0 

ELSEIF  (Out  (I,  J).GT. 65535)  THEN 
Out  (I,  J)  =  65535 
END  IF 

100  CONTINUE 

RETURN 

END 
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Y. _ Temporal  Filtering 

A.  Description 

The  temporal  filtering  algorithm  provides  a  pixel-by-pixel  infinite  impulse  response  (HR) 
filtering  of  the  sequence  of  images  which  are  sent  through  the  process.  This  implementation 
allows  up  to  a  fourth-order  filter. 

B.  Data  Fields 


Action 

Field  Header 

Field  Data 

Data  Type 

input 

Dimensions 

FPA  dimension 

Integer 

E&39HIH 

Temporal  Filtering  Limits 

Lower  limits 

Integer 

Integer 

modify 

Pixel  Data 

Pixel  data  array 

Integer 

[1.. Dimension] 

[1..  Dimension] 

C.  Module  Listing 

PROGRAM  TEMPORALFILTER 
C 

C  Temporal  Filter  Test  Module 

C 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 

C  400  Tenth  St.  CRB  390 

C  Atlanta,  GA  30332-0540 

C  Contact:  Andrew  Henshaw  (404)894-2521 

C  . 

C  conforms  to  the  Ga.  Tech  Algorithm  Evaluation  Data  Format 
C 

C  Fortran  translation  of  Occam  code 

C  Steve  Gieseking 

C  Roy  Melton 

C 

C  Harness  written  by  A.  M.  Henshaw  Jan  23,  1990 

C  Using  Microsoft  Fortran 

C 

CHARACTER* (* )  Dim,  Pixels,  Limits 
C 

C  Valid  Section  Headers 

PARAMETER  (Dim  =' Dimensions ' ) 

PARAMETER  (Pixels* ' Pixel  Data') 

PARAMETER  (Limits* ' Temporal  Filtering  Limits') 

C 

PARAMETER  (maxSize=64)  !  maximum  FPA  size 
PARAMETER  (TF_CONTROL_SIZE  =  24) 

INTEGER  n,  lower,  upper 

INTEGER  in (maxSize,maxSize) ,  out (maxSize,  maxSize) 

INTEGER  X  (maxSize,  maxSize,  2,  2) 

INTEGER  C  ( TF_C0NTR0L_S I ZE ) 

CHARACTER  header*72 
LOGICAL  runFlag 
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DATA  lower  /0/ 

DATA  upper  / 327 67/ 

WRITE  (6,*)  '%  Processed  by  Temporal  Filtering  Module.' 

runFlag  -  .TRUE. 

DO  WHILE  (runFlag) 

READ  (5,1000)  header 
1000  FORMAT  (A72)  ' 

IF  (header. EQ. Dim)  THEN 
READ  (5,*)  n 
WRITE  (6,*)  Dim 
WRITE  (6,  *)  n 

ELSE  IF  (header. EQ. Limits)  THEN 
READ  (5,*)  lower,  upper 
WRITE  (6,*)  Limits 
WRITE  (6,*)  lower,  upper 
ELSE  IF  (header. EQ. Pixels)  THEN 

READ  (5,*)  ( (in (row, col) , col=l, n) , row=l,  n) 

CALL  CalculateFilterControls  (C,  Lower,  Upper) 

CALL  TempFilt  (In,  Out,  X,  C,  N) 

WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (out (row, col) ,col=l,n) , row=l,n) 

ELSE  IF  (header . EQ . 1  End 1 )  THEN 
WRITE  (6,*)  'End' 
runFlag  =  .FALSE. 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 


END 


C***Filter  Control  calculation*************************************** 
SUBROUTINE  CalculateFilterControls  (Control,  Lower,  Upper) 


PARAMETER  (TF_CONTROL_SIZE  -  24) 

INTEGER  Control  (TF_CONTROL_SIZE) 
INTEGER  Lower,  Upper 
INTEGER  I,  J 

DO  110  1=0,  12,  12 
DO  100  J  -  1,  8 

Control  (J  +  I)  =1 
100  CONTINUE 


Control 
Control 
Control 
Control 
110  CONTINUE 


(9  +  1)  =3 
(10  +  I)  =  1 

(11  +  I)  =  Upper 

(12  +  I)  =  Lower 


RETURN 

END 


C***Temporal  Filter************************************************** 
SUBROUTINE  TempFilt  (In,  Out,  X,  C,  N) 


PARAMETER  (MAXSIZE  =  64) 
PARAMETER  (TF_A0  =  1) 
PARAMETER  (TF  A1  =  2) 
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PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 


(TF_A2  =  3) 

(TF_B0  =  4) 

(TF_B1  -  5) 

<TF_B2  =  6) 

( TF_SCALE_ST ATE  =  7) 
(TF_SCALE_OUTPUT  =  8) 
(TF_UPPER_LIMIT_STATE  =  9) 

( T  F_LOWE  R_L IM I T_S TATE  =  10) 
(TF_UPPER_LIMIT_OUTPUT  =  11) 
(TF_LOWER_LIMIT_OUTPUT  =  12) 
(TF_CONTROL  SIZE  =24) 


INTEGER  In  (MAXSIZE,  MAXSIZE),  Out  (MAXSIZE,  MAXSIZE) 
INTEGER  X  (MAXSIZE,  MAXSIZE,  2,  2) 

INTEGER  C  (TF_CONTROL_SIZE) 

INTEGER  N 

INTEGER  I,  J,  K,  L,  Ptr,  Value,  XNew,  YNew 

DO  10  I  =  1,  MAXSIZE 
DO  10  J  =  1,  MAXSIZE 
DO  10  K  -  1,  2 
DO  10  L  =  1,  Z 

X  (I,  J,  K,  L)  =  0 

10  CONTINUE 

DO  30  I  =  1,  N 
DO  30  J  =  1,  N 
Value  =  In  (I,  J) 

DO  20  K  -  1,2 

Ptr  =  (K  -  1)  *  12 


XNew  =  ( (C  (Ptr 

+  TF  A0) 

*  Value) 

+ 

+ 

(C  (Ptr 

+  TF  Al) 

*  X  (I, 

J , 

K, 

D)  + 

+ 

(C  (Ptr 

+  TF  A2) 

*  X  (I, 

J, 

K, 

2)  )  ) 

/ 

+ 

C  (Ptr  + 

TF  SCALE 

STATE) 

YNew  =  ( (C  (Ptr 

+  TF  B0)' 

'*  Value) 

+ 

+ 

(C  (Ptr 

+  TF  Bl) 

*  X  (I, 

J, 

K, 

1)  )  + 

+ 

(C  (Ptr 

+  TF  B2) 

*  X  (I, 

J, 

K, 

2)  )  ) 

/ 

+ 

C  (Ptr  + 

TF  SCALE 

OUTPUT) 

X  (I,  J,  K,  2)  =  X  (I,  J,  K,  1) 

IF  (XNew.GT.C  (Ptr  +  TF_UPPER_LIMIT_STATE) )  THEN 
X  (I,  J,  K,  1)  =  C  (Ptr  +  TF_UPPER_LIMIT_STATE) 
ELSEIF  (XNew.LT.C  (Ptr  +  TF_LOWER_LIMIT_STATE) )  THEN 
X  (I,  J,  K,  1)  =  C  (Ptr  +  TF_LOWER_L IMI T_S  TATE ) 

ELSE 

X  (I,  J,  K,  1)  =  XNew 
END  IF 

IF  (YNew.GT.C  (Ptr  +  TF_UPPER_LIMIT_OUTPUT) )  THEN 
Value  =  C  (Ptr  +  TF_UPPER_LIMIT_OUTPUT) 

ELSEIF  (YNew.LT.C  (Ptr  +  TF_LOWER_LIMIT_OUTPUT) )  THEN 
value  =  C  (Ptr  +  TF_LOWER_LIMIT_OUTPUT) 

ELSE 

Value  =  YNew 
END  IF 

20  CONTINUE 

Out  (I,  J)  =  Value 
30  CONTINUE 


RETURN 

END 
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YL _ Thresholding 

The  thresholding  algorithm  is  used  to  partition  the  image  into  points  which  are  of  interest  and 
those  that  are  not  of  interest  Pixels  are  zeroed  if  they  are  not  of  interest.  A  pixel  is  passed  if  the 
intensity  is  above  a  calculated  lower  threshold  value  and  below  a  fixed  upper  threshold  value. 
The  lower  threshold  supports  two  of  the  modes  which  are  in  the  Georgia  Tech  VLSI  design. 
This  includes  a  simple,  fixed  threshold  and  an  adaptive  threshold  based  on  the  average  and  first 
central  absolute  moment  of  the  surrounding  eight  pixels. 

YIL _ Simple  Thresholding 

1.  Data  Fields 


Action 

Field  Header 

Field  Data 

Data  Type 

input 

Dimensions 

FPA  dimension 

‘Integer 

BSsSSSHHI 

Simple  Thresholding  Limits 

Lower  limit 

Integer 

Integer 

:E53I^3HHI 

Simple  Thresholding  Statistics 

Lower  limit  used 

Integer 

Integer 

Count  of  pixels 
exceeding  limit 

Integer 

modify 

Pixel  Data 

Pixel  data  array 

Integer 

[1.. Dimension] 
ri..  Dimension] 

2.  Module  Listing 

PROGRAM  STHRESH 
C 

C  Simple  Thresholding  Test  Module 

C 

C  conforms  to  the  Ga .  Tech  Algorithm  Evaluation  Data  Format 
C 

C  Fortran  translation  of  Occam  code 

C  Steve  Gieseking 

C  Andrew  Henshaw 

C 

C  Harness  written  by  Andrew  Henshaw  Jan  23,  1990 

C  Using  Microsoft  Fortran 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 

C 

CHARACTER* <*)  Dim,  Pixels,  Limits 
C 

C  Valid  Section  Headers 

C 

PARAMETER  (Dim  = 'Dimensions ' ) 

PARAMETER  (Pixels= ' Pixel  Data') 

PARAMETER  (Limits= ' Simple  Thresholding  Limits') 
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PARAMETER  (maxSize=128)  !  maximum  FPA  size 

INTEGER  n,  count,  lower,  upper 

INTEGER  in (maxSize, maxSize) ,  out(maxSize,  maxSize) 

CHARACTER  header*72 
LOGICAL  runFlag 

DATA  lower  /0/  !  default  values 

DATA  upper  / 327 67/ 

WRITE  (6,*)  '%  Processed  by  Simple  Thresholding  module.' 
runFlag  «  . TRUE . 

DO  WHILE  (runFlag) 

READ  (5,1000)  header 
1000  FORMAT  (A72) 

IF  (header .EQ. Dim)  THEN 
READ  (5,*)  n 
WRITE  (6,*)  Dim 
WRITE  (6,*)  n 

ELSE  IF  (header. EQ. Limits)  THEN 
READ  (5,*)  lower,  upper 
WRITE  (6,*)  Limits 
WRITE  (6,*)  lower,  upper 
ELSE  IF  (header. EQ. Pixels)  THEN 

READ  (5,*)  ( (in (row, col) , col=l, n) , row=l, n) 

CALL  SmpThrsh  (n,  lower,  upper,  count,  in,  out) 

WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (out (row, col) ,col=l,n) , row=l,  n) 

ELSE  IF  (header. EQ. 'End' )  THEN 
WRITE  (6,*)  'End' 
runFlag  -  .FALSE. 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 

END 

c* ******************************************************************* 

SUBROUTINE  SmpThrsh  (n,  lower,  upper,  count,  in,  out) 

PARAMETER  (maxSize=*64) 

INTEGER  n,  lower,  upper,  count 
INTEGER  in (maxSize,  maxSize) 

INTEGER  out (maxSize,  maxSize) 

INTEGER  row,  col,  pixel 


count  =  0 
DO  30  row  =  1,  n 
DO  30  col  =  1,  n 
pixel  =  in (row, col) 

IF  ( (pixel .GE . lower) .AND . (pixel . LE . upper) )  THEN 
count  =  count  +  1 
out (row, col)  =  pixel 
ELSE 

out (row, col)  =  0 
END  IF 

30  CONTINUE 

C  Put  Statistics  onto  data  stream 
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WRITE  (6,*)  'Simple  Thr4sholding  Statistics' 

WRITE  (6,*)  lower,  upper,  count 

RETURN 

END 
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Vm.  Adaptive  Thresholding 
1.  Data  Fields 


Action 

Field  Header 

Field  Data 

Data  Type 

input 

Dimensions 

FPA  dimension 

Integer 

E339HHI 

Adaptive  Thresholding  Parameters 

Integer 

kl 

Integer 

k2 

Integer 

k3 

Integer 

Scale 

Integer 

Adaptive  Thresholding  Statistics 

Integer 

Count  of  pixels 
exceeding  limit 

Integer 

modify 

Pixel  Data 

Pixel  data  array 

Integer 

[1.. Dimension] 

[L.  Dimension] 

2.  Module  Listing 

PROGRAM  ADTHRESH 

C 

C  Adaptive  Thresholding  Test  Module 

C 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 

C  400  Tenth  St.  CRB  390 

C  Atlanta,  GA  30332-0540 

C  Contact:  Andrew  Henshaw  (404)894-2521 

C 

C  conforms  to  the  Ga.  Tech  Algorithm  Evaluation  Data  Format 
C 

C  Fortran  translation  of  Occam  code 

C  Steve  Gieseking 

C  Roy  Melton 

C  Harness  written  by  A.  M.  Henshaw  Jan  23,  1990 

C  Using  Microsoft  Fortran 

C 

CHARACTER* (*)  Dim,  Pixels 
PARAMETER  (Dim  “'Dimensions  ') 

PARAMETER  (Pixels= ' Pixel  Data  ') 

PARAMETER  (Parms= ' Adaptive  Thresholding  Parameters') 
PARAMETER  (maxSize=64) 

INTEGER  n,  count,  kl,  k2,  k3,  scale,  sum,  upper 
INTEGER  in (maxSize, maxSize) ,  out(maxSize,  maxSize) 
CHARACTER  header*72 
LOGICAL  runFlag 
DATA  kl  /!/ 

DATA  k2  /0/ 

DATA  k3  /0/ 

DATA  scale  / 8 / 


-18- 


Adaptive  Thresholding 


DATA  upper  /327 67/ 

WRITE  <6,*)  '%  Processed  by  Adaptive  Thresholding  module.' 
runFlag  =  .TRUE. 

DO  WHILE  (runFlag) 

READ  (5,1000)  header 
1000  FORMAT  (A72) 

IF  (header. EQ. Dim)  THEN 
READ  (5,*)  n 
WRITE  (6,*)  Dim 
WRITE  (6,*)  n 

ELSE  IF  (header. EQ.Parms)  THEN 

READ  (5,*)  upper,  kl,  k2,  k3,  scale 
WRITE  (6,*)  Parms 

WRITE  (6,*)  upper,  kl,  k2,  k3,  scale 
ELSE  IF  (header. EQ. Pixels)  THEN 

READ  (5,*)  ( (in (row, col) , col-1, n) , row-1, n) 

CALL  AdThrsh  (n,  upper,  count,  sum, 

+  kl,  k2,  k3,  scale,  in,  out) 

WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (out (row, col) , col-1, n) , row=l, n) 

ELSE  IF  (header. EQ. 'End' )  THEN 
WRITE  (6,*)  'End' 
runFlag  =  . FALSE . 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 

END 

C* ******************************************************************* 

SUBROUTINE  AdThrsh  (N,  Upper,  Count,  Sum,  Kl,  K2,  K3,  Scale, 

+  In,  Out) 

PARAMETER  (maxSize=64) 

INTEGER  N,  Upper,  Count,  Sum,  Kl,  K2,  K3,  Scale 
INTEGER  In(maxSize,  maxSize) 

INTEGER  Out (maxSize,  maxSize) 

INTEGER  Average,  I,  J,  K,  L,  Lower,  Stat 


Count  =  0 
Sum  =  0 

DO  30  I  =  1,  N 
DO  30  J  =  1,  N 

IF  ( ( (I.EQ.l) .OR. (I.EQ.N) ) .OR. ( (J.EQ.l) .OR. (J.EQ.N) ) )  THEN 
Out  (I,  J)  =0 
ELSE 

Average  =  In  (1-1,  J-l)  +  In  (1-1,  J)  +  In  (1-1,  J+l)  + 

+  In  (I,  J-l)  +  In  (I,  J+l)  + 

+  In  (1+1,  J-l)  +  In  (I+l,  J)  +  In  { 1+1 ,  J+l) 

Stat  =  0 

DO  10  K - 1,  1 

DO  10  L - 1,  1 

IF  ( (K.NE.0) .AND. (L.NE.0) )  THEN 

Stat  =  Stat  +  ABS  ((In  (I+K,  J+L)  *  8)  -  Average) 

END  IF 

10  CONTINUE 

Lower  =  ((Average  *  Kl)  +  (Stat  *  K2)  +  K3)  /  Scale 
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Sum  =  Sum  +  Stat 

IF  ((In  (I,  J) .GE. Lower) .AND. (In  (I,  J) . LE . Upper ) )  THEN 
Out  (I,  J)  =  In  (I,  J) 

Count  =  Count  +  1 
ELSE 

Out  (I,  J)  =  0 
ENDIF 
ENDIF 

30  CONTINUE 

C  Put  Statistics  onto  data  stream 

WRITE  (6,*)  'Adaptive  Thresholding  Statistics' 

WRITE  (6,*)  Upper,  Count 

RETURN 

END 
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EL _ Clustering  &  Centroiding 

A.  Description 

The  clustering  algorithm  forms  connected  sets  of  pixels  based  on  the  surrounding  pixels.  Two 
pixels  are  elements  of  the  same  cluster  of  pixels  if  they  are  one  of  the  eight  nearest  neighbors  of 
each  other.  The  centroiding  algorithm  calculates  the  area  centroid  and  the  intensity  weighted 
centroid  of  the  clusters  specified  by  the  clustering  algorithm. 

B.  Data  Fields 


Action 

Field  Header 

Field  Data 

Data  Type 

input 

Dimensions 

FPA  dimension 

Integer 

input 

Pixel  Data 

Pixel  data  array 

Integer 

[1.. Dimension] 

[1..  Dimension] 

Clusters 

Ouster  count 

Integer 

output 

Centroids 

Vector  of  the 
following  repeated 
Ouster  count  times 

Integer 

Integer 

Intensity  centroid  (X) 

Integer 

Intensity  centroid  (Y) 

Integer 

Area  in  pixels 

Integer 

Total  cluster  intensity 

Integer 

C.  Module  Listing 

PROGRAM  CENTROID 
C 

C  Clustering  and  Centroiding  Test  Module 

C 

C  Computer  Engineering  Research  Laboratory 

C  Georgia  Institute  of  Technology 

C  400  Tenth  St.  CRB  390 

C  Atlanta,  GA  30332-0540 

C  Contact:  Andrew  Henshaw  (404)894-2521 

C 

C  conforms  to  the  Ga .  Tech  Algorithm  Evaluation  Data  Format 

C 

C  Fortran  translation  of  Occam  code 

C  Steve  Gieseking 

C  Roy  W.  Melton  Feb  12,  1990 
C 

C  Harness  written  by  Andrew  Henshaw  Jan  23,  1990 

C  Using  Microsoft  Fortran 

C 

CHARACTER* (*)  Dim,  Pixels 
PARAMETER  (Dim  =' Dimensions  ') 

PARAMETER  (Pixels= ' Pixel  Data  ') 

PARAMETER  (MAX_SIZE=64) 

PARAMETER  (MAX_CLUSTERS=1024) 
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INTEGER  ClusterCount,  N 

INTEGER  Frame  (MAX_SIZE,  MAX_SIZE) 

INTEGER  Clusters  (MAX_CLUSTERS,  6) 

CHARACTER  header*72 
LOGICAL  runFlag 

WRITE  (6,*)  '%  Processed  by  Centroid  Image  Module.'  • 
runFlag  =  .TRUE. 

DO  WHILE  (runFlag) 

READ  (5,1000)  header 
1000  FORMAT  (A72) 

IF  (header. EQ. Dim)  THEN 
READ  (5,*)  N 
WRITE  (6,*)  Dim 
WRITE  (6,*)  N 

ELSE  IF  (header. EQ. Pixels)  THEN 

READ  (5,*)  ( (Frame (row, col) , col-1, N) , row-1, N) 

CALL  Centroidlmage  (Frame,  Clusters,  N,  ClusterCount) 

WRITE  (6,*)  Pixels 

WRITE  (6,*)  ( (Frame (row, col) , col-1,  N) ,  row-1, N) 

WRITE  (6,*)  'Clusters' 

WRITE  (6,*)  ClusterCount 
IF  (ClusterCount. GT.0)  THEN 
WRITE  (6,*)  'Centroids' 

WRITE  (6,*)  ((Clusters  (row,  col),  col-1,6), 

+  row-1, ClusterCount) 

ENDIF 

ELSE  IF  (header. EQ. 'End' )  THEN 
WRITE  (6,*)  'End' 
runFlag  -  . FALSE . 

ELSE 

WRITE  (6,*)  header 
END  IF 
END  DO 

END 


SUBROUTINE  Centroidlmage  (Frame,  CData,  N,  ClusterCount) 


PARAMETER  (MAX_SIZE«64 ) 

PARAMETER  (MAX_CLUSTERS=102 4 ) 

PARAMETER  (CSum  -  1) 

PARAMETER  (CSumX  -  2) 

PARAMETER  (CSumY  =  3) 

PARAMETER  (ISum  =  4) 

PARAMETER  (ISumX  =  5) 

PARAMETER  (ISumY  =  6) 

PARAMETER  (ACoorX  =  1) 

PARAMETER  (ACoorY  =  2) 

PARAMETER  (ICoorX  =  3) 

PARAMETER  (ICoorY  =  4) 

PARAMETER  (Area  =  5) 

PARAMETER  (Intensity  =  6) 

INTEGER  Frame  (MAX_SIZE,  MAX_SIZE) 

INTEGER  CData  (MAX_CLUSTERS,  6) 

INTEGER  N,  ClusterCount 

INTEGER  CO,  Cl,  CNM1,  CN,  CNP1,  FinalCluster,  I,  J 
INTEGER  Cluster  (MAX_SIZE  +  1),  Temp  (6) 

INTEGER  Reassign 
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DIMENSION  Reassign  (0 :MAX_CLUSTERS-1) 

DO  10  1=1, N+l 
Cluster  (I)  =0 
10  CONTINUE 

FinalCluster  =  0 
Reassign  (0)  =  0 

DO  50  1=1, N 

C  /*  Initialize  Row  */ 

Cl  =  0 

CNP1  =  0 

CN  =  Cluster  (1) 

20  IF  (CN.NE. Reassign  (CN) )  THEN 
CN  =  Reassign  (CN) 

GOTO  20 
END  IF 

DO  40  J  =  1,  N 

CNM1  =  Cluster  (J+l) 

30  IF  (CNM1.NE. Reassign  (CNM1) )  THEN 
CNM1  =  Reassign  (CNM1) 

GOTO  30 
END  IF 

IF  (Frame  (I,  J).EQ.0)  THEN 
CO  =  0 

ELSEIF  (C1.NE.0)  THEN 
C  /*  Add  Pixel  to  Cluster  */ 

IF  ( (CNM1.NE.0) .AND. (CNM1.NE.C1) )  THEN 
C  /*  Merge  Cl,  CNM1  */ 

CData  (Cl,  CSum)  =  CData  (Cl,  CSum)  + 

+  CData  (CNM1,  CSum)  +  1 

CData  (Cl,  CSumX)  =  CData  (Cl,  CSumX)  + 

+  CData  (CNM1,  CSumX)  +  J 

CData  (Cl,  CSumY)  =  CData  (Cl,  CSumY)  + 

+  CData  (CNM1,  CSumY)  +  I 

CData  (Cl,  ISum)  =  CData  (Cl,  ISum)  + 

+  CData  (CNM1,  ISum)  +  Frame  (I,  J) 

CData  (Cl,  ISumX)  =  CData  (Cl,  ISumX)  + 

+  CData  (CNM1,  ISumX)  + 

+  (Frame  (I,  J)  *  J) 

CData  (Cl,  ISumY)  =  CData  (Cl,  ISumY)  + 

+  CData  (CNM1,  ISumY)  + 

+  (Frame  (I,  J)  *  I) 

CData  (CNM1,  CSum)  =  0 
Reassign  (CNM1)  =  Cl 
CNM1  =  Cl 
CN  =  Reassign  (CN) 

CO  =  Cl 

ELSE 

/*  Add  to  Cl  */ 

CData  (Cl,  CSum)  =  CData  (Cl,  CSum)  +  1 
CData  (Cl,  CSumX)  =  CData  (Cl,  CSumX)  +  J 
CData  (Cl,  CSumY)  =  CData  (Cl,  CSumY)  +  I 
CData  (Cl,  ISum)  =  CData  (Cl,  ISum)  +  Frame  (I,  J) 
CData  (Cl,  ISumX)  =  CData  (Cl,  ISumX)  + 

(Frame  (I,  J)  *  J) 

CData  (Cl,  ISumY)  =  CData  (Cl,  ISumY)  + 

(Frame  (I,  J)  *  I) 

CO  =  Cl 
END  IF 

ELSEIF  (CN.NE. 0)  THEN 
C  /*  Add  to  CN  */ 
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CData  (CN, 
CData  (CN, 
CData  (CN, 
CData  (CN, 
CData  (CN, 
CData  (CN, 
CO  =  CN 


CSum)  =  CData  (CN,  CSum)  +  1 

CSumX)  =  CData  (CN,  CSumX)  +  J 

CSumY)  =  CData  (CN,  CSumY)  +  I 

ISum)  =  CData  (CN,  ISum)  +  Frame  (I,  J) 

ISumX)  =  CData  (CN,  ISumX)  +  (Frame  (I, 

ISumY)  =  CData  (CN,  ISumY)  +  (Frame  (I, 


J) 

J) 


J) 

I) 


C 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


ELSE IF  (CNM1.NE.0)  THEN 

IF  ( (CNP1.NE.0) .AND. (CNP1.NE.CNM1) )  THEN 
/*  Merge  CNM1,  CNP1  */ 

CData  (CNM1,  CSum)  =  CData  (CNM1,  CSum)  + 

CData  ( CNP 1 ,  CSum)  +  1 
CData  (CNM1,  CSumX)  =  CData  (CNM1,  CSumX)  + 

CData  ( CNP 1 ,  CSumX)  +  J 
CData  (CNM1,  CSumY)  =  CData  (CNM1,  CSumY)  + 

CData  (CNP1,  CSumY)  +  I 
CData  (CNM1,  ISum)  =  CData  (CNM1,  ISum)  + 

CData  (CNP1,  ISum)  +  Frame  (I,  J) 
CData  (CNM1,  ISumX)  «=  CData  (CNM1,  ISumX)  + 

CData  ( CNP 1 ,  ISumX)  + 

(Frame  (I,  J)  *  J) 

CData  (CNM1,  ISumY)  =  CData  (CNM1,  ISumY)  + 

CData  (CNP1,  ISumY)  + 

(Frame  (I,  J)  *  I) 

CData  (CNP1,  CSum)  =  0 
Reassign  (CNP1)  =  CNM1 
CO  =  CNM1 


C 


+ 

+ 


ELSE 

/*  Add  to  CNM1  */ 

CData  (CNM1 ,  CSum)  -  CData  (CNM1,  CSum)  + 
CData  (CNM1 ,  CSumX)  -  CData  (CNM1,  CSumX) 
CData  (CNM1,  CSumY)  =  CData  (CNM1,  CSumY) 
CData  (CNM1,  ISum)  =  CData  (CNM1,  ISum)  + 
CData  (CNM1,  ISumX)  =  CData  (CNM1,  ISumX) 

(Frame  (I,  J)  *  J) 

CData  (CNM1,  ISumY)  =  CData  (CNM1,  ISumY) 

(Frame  (I,  J)  *  I) 


CO  = 
END  IF 


CNM1 


1 

+  J 
+  I 

Frame  (I,  J) 
+ 


ELSEIF  (CNP1.NE.0)  THEN 
C  /*  Add  to  CNP1  */ 

CData  (CNP1,  CSum)  =  CData  (CNP1,  CSum)  +  1 

CData  (CNP1,  CSumX)  =  CData  (CNP1,  CSumX)  +  J 

CData  ( CNP 1 ,  CSumY)  =  CData  (CNP1,  CSumY)  +  I 

CData  ( CNP 1 ,  ISum)  =  CData  (CNP1,  ISum)  +  Frame  (I,  J) 

CData  (CNP 1 ,  ISumX)  =  CData  (CNP1,  ISumX)  + 

+  (Frame  (I,  J)  *  J) 

CData  (CNP1,  ISumY)  =  CData  (CNP1,  ISumY)  + 

+  (Frame  (I,  J)  *  I) 

CO  =  CNP1 


C 


ELSE 

/*  New  Cluster  */ 

FinalCluster  =  FinalCluster  +  1 


CO  =  FinalCluster 


CData 

CData 

CData 

CData 

CData 

CData 


(CO,  CSum)  =  1 

(CO,  CSumX)  =  J 

(CO,  CSumY)  =  I 

(CO,  ISum)  =  Frame  (I,  J) 

(CO,  ISumX)  =  Frame  (I,  J) 

(CO,  ISumY)  =  Frame  (I,  J) 


*  J 

*  I 
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Reassign  (CO)  =  CO 
END  IF 


Cluster  (J)  -  CO 


C  /*  Update  for  next  column  */ 

Cl  =  CO 
CNP1  -  CN 
CN  =  CNM1 


40  CONTINUE 
50  CONTINUE 


/*  Output  Centroids  */ 
ClusterCount  =  0 
DO  70  I=l,FinalCluster 

IF  (CData  (I,  CSum).NE.O)  THEN 
/*  Valid  Cluster  */ 


+ 

Temp 

(ACoorX)  -  (CData  (I,  CSumX)  + 
ISHFT  (CData  (I,  CSum) , 

-1) 

) 

/ 

CData 

(I,  CSum) 

+ 

Temp 

(ACoorY)  -  (CData  (I,  CSumY)  + 
ISHFT  (CData  (I,  CSum) , 

-1) 

) 

/ 

CData 

(I,  CSum) 

+ 

Temp 

(ICoorX)  =  (CData  (I,  ISumX)  + 
ISHFT  (CData  (I,  ISum) , 

.-1) 

) 

/ 

CData 

(I,  ISum) 

+ 

Temp 

(ICoorY)  -  (CData  (I,  ISumY)  + 
ISHFT  (CData  (I,  ISum), 

-1) 

) 

/ 

CData 

(I,  ISum) 

Temp 

Temp 

(Area)  -  CData  (I,  CSum) 
(Intensity)  -  CData  (I,  ISum) 

60 


ClusterCount  -  ClusterCount  +  1 
DO  60  J=l,6 

CData  (ClusterCount,  J)  =  Temp 
CONTINUE 
END  IF 


(J) 


70  CONTINUE 


RETURN 

END 


-25- 


